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. Abstract. We perform the Bayesian inference of a GARCH model by the Metropolis-Hastings 

algorithm with an adaptive proposal density. The adaptive proposal density is assumed to be 
the Student's t-distribution and the distribution parameters are evaluated by using the data 
' sampled during the simulation. We apply the method for the QGARCH model which is one 

of asymmetric GARCH models and make empirical studies for for Nikkei 225, DAX and Hang 
indexes. We find that autocorrelation times from our method are very small, thus the method 
^ ' is very efficient for generating uncorrelated Monte Carlo data. The results from the QGARCH 

' [~l I model show that all the three indexes show the leverage effect, i.e. the volatility is high after 

^ |i negative observations. 



> , 

^ . 1. Introduction 

In empirical finance volatility of asset returns is an important value to measure risk. In order 
Ps| ■ to forecast future volatility it is desirable to use appropriate models which have the properties 

QQ , of volatility of asset returns. Many empirical studies suggest that the distribution of asset 

I returns is leptokurtic. Furthermore the volatility is not constant, but changes over time. There 

' are periods when volatility is very high or very low. This property of the volatility is called 

volatility clustering. 

, ^ ' The Autoregressive Conditional Heteroscedasticity (ARCH) model[l] and its generalization, 

^ . the Generalized ARCH (GARCH) model[2] are designed to capture the property of the volatility 

^ I clustering. Moreover the distributions of returns from those models show fat-tailed distributions 

and they are suggested to be Student's t-(Tsallis) distributions [3lllj. There are many extensions 
of the GARCH model to include additional properties of the volatility. An example of the 
properties of the volatility is that the volatility response is high after negative news (returns), 
which is known as the leverage effect, first observed by Black In order to cope with this 
fact, some modelsjU EJ El El HOl [H] which introduce asymmetry into the volatility response 
function are proposed. In this study among them we focus on the Quadratic GARCH(QGARCH) 
model ^lOt illj which adds an additional term proportional to a return to the volatility response 
function. 

To utilize the GARCH models we need to infer GARCH parameters from financial time series 
data. In general the Maximum Likelihood (ML) estimation is favored to the inference of GARCH 
models. Although implementation of the ML method is straightforward, there exist practical 
difficulties in estimating GARCH parameters by the ML technique. The model parameters must 
be positive to ensure a positive volatility and the stationarity condition for volatility is also 



required. The ML method with such requirements is performed via a constrained optimization 
technique which can be cumbersome. Forethermore the output of the ML method is often 
sensitive to starting values. 

Another estimation technique is the Bayesian inference which does not have the difficulties 
seen in the ML method. Usually the Bayesian inference is performed by Markov Chain 
Monte Carlo (MCMC) methods which have been common in the recent computer development. 
Various MCMC methods for the Bayesian inference of the GARCH models have been 
proposeddl [HEKISKISKITKISKIDj. In a survey on the MCMC methods of the GARCH 
models[T7] it is shown that Acceptance- Rejection/ Metropolis-Hastings (AR/MH) algorithm 
with a multivariate Student's t-distribution works better than the other algorithms. The 
multivariate Student's t-distribution is used as a proposal density of the MH algorithm and the 
parameters to specify the Student's t-distribution are determined by the Maximum Likelihood 
(ML) technique. Recently an alternative method to estimate those parameters without relying 
on the ML technique was proposed [201 Ell ES]- The method is called "adaptive construction 
scheme", where the parameters of the multivariate Student's t-distribution are determined by 
using the pre-sampled data by an MCMC method. And the parameters are updated adaptively 
during the MCMC simulation. 

The adaptive construction scheme was tested for GARCH and QGARCH models [20 1 12 H l22] 
and it is shown that the adaptive construction scheme can significantly reduce the correlation 
between sampled data. In this paper first we describe the adaptive construction scheme for the 
GARCH models. Then we make empirical studies with the QGARCH model for three major 
stock indexes, Nikkei 225, DAX and Hang Seng. 

2. GARCH Model 

Let xt be an asset return observed at time t. We transform xt to yt as 

yt = xt- X, (1) 

where x is the average over observations, i.e. x 
assumed to be decomposed as 

yt = crtet, 

where et is an identically distributed random variable with zero mean and unit variance. The 
distribution of et is usually assumed to be a normal (Gaussian) or a leptokurtic one. In this 
study we assume that the distribution is a normal one, i.e. ~ A^(0, 1). In the original GARCH 
model the volatility at is assumed to change over time as 

q p 

<^t =^ + Yl ^iyt-i + f^i^t-i > (3) 

i=l i=l 

and more specifically this model is stated as the GARCH(p,q) model. Since in empirical studies 
with the AIC analysis small numbers are often chosen for p and q we focus on the GARCH(1,1) 
model in this study and write it as the GARCH model. Now eq.® is written as 

= IV + ay^_^ + pa^_^. (4) 

The volatility response function of eq.dH) is symmetric under positive or negative observations 
yt- However some asset returns are known to show the leverage effect, i.e. the volatility is higher 
after negative observations than after positive ones. Several extended GARCH models have been 



= TT^^^j. In the GARCH model yt is 

i=i 

(2) 



proposed to introduce asymmetry into the volatility response function. In this study we use the 
QGARCH model [ini ttl] given by 

af = u + jyt-i + ayi_i + /?cri_i , (5) 



where the additional term, jyt^i introduces the asymmetry into the model. The QGARCH 
model includes 4 parameters (cj, a, f3, 7) which have to be determined from financial data. 



3. Maximum Likelihood Estimation 

The likelihood function of the GARCH models is written as 

Liy\e) = nr=i^^exp (-^), (6) 

where y stands for the time series data of n observations, y = {yi, ■■■■yn) and 9 stands for GARCH 
parameters, e.g. for the QGARCH model 9 = (u;,a,/3,7). 

In the ML estimation the parameters are determined by maximizing the log likelihood 
lnL{y\9), 

-.71 n 2 

In L{y\9) = --^ln(2vra?) - E (7) 

7 i i 



4. Bayesian Inference 

From the Bayes' theorem we obtain the posterior density vr(0|?/) which is a probability 
distribution of parameters 9 as 

n{9\y)^L{y\9)7T{9), (8) 

where Tr{9) is the prior density for 9. Since we do not know the functional form of vr(0) we make 
an assumption for it. Here we assume that it{9) is constant. Using TT{9\y), 9 is estimated by 

(^) = I / 07T{9\y)d9, (9) 

where 

Z = J 7r{9\y)d9. (10) 

In general eq.(l9]) can not be performed analytically. Instead of performing the integral of eq.® 
we evaluate eq.Q by the MCMC technique described in the next section. In the MCMC 
evaluation the normalization constant Z will be irrelevant. 



5. Markov Chain Monte Carlo 

The Metropolis algorithm is an MCMC technique first introduced by Metropolis et a/.[23j, 
which is designed to estimate an integral such as eq.([9]) numerically. The Metropolis-Hastings 
(MH) algorithm |24j is a generalization of the original Metropolis algorithm. Let us consider 
to evaluate eq.© by the MH algorithm. In general the functional form of 7r{9\y) may be too 
complicated to generate random variables according to TT(9\y). In the MH algorithm we draw 
random variables from a proposal density which is simple enough to generate random variables. 
The basic procedure of the MH algorithm is given as follows. 

(1) First we set an initial value and i = 1. 

(2) Then we generate a new value 9i from a certain proposal density g{9i\9i-i). 



(3) We accept the candidate 6i with a probabihty of PMH(di^i^6i) where 



Pmh{ 



mm 



1, 



When 9i is rejected we keep Oi-i, i.e. 9i = Oi^i. 

(4) Go back to (2) with an increment of i = i + 1. 

When the proposal density does not depend on the previous value, i.e. g{0i\9i-i) = g{ 
obtain 



PMH{9i~i,0i) = min 



1, 



(11) 



we 



(12) 



P{Oi-i) giOi) 

For a symmetric proposal density gi6i\6i^i) = g{6i-i\9i)^ eq. (fTT]) reduces to the Metropolis 
algorithm and the Metropolis accept probability is given by 

-PMetro (^i-l 5 6'j) = min 1, ,f ^\ . 

Eq.([9]) is evaluated as an average over the data sampled by the MCMC algorithm. 

hm 
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i=l 



where N is the number of the sampled data. For N independent data the statistical error is 
proportional to However the data sampled by MCMC methods are not independent. For 

correlated data, the statistical error is proportional to where r is the autocorrelation time 
between the sampled data. Thus in order to have a small statistical error without increasing the 
number of sampled data, it is important to take an MCMC method which generates uncorrelated 
data. 



6. Adaptive Construction Scheme 

The efficiency of the MH algorithm depends on how we choose the proposal density. By choosing 
an adequate proposal density for the MH algorithm one can reduce the correlation between the 
sampled data. The posterior density of GARCH parameters often resembles to a Gaussian-like 
shape. Thus one may choose a density similar to a Gaussian distribution as the proposal density. 
Such attempts have been done by Mitsui, Watanabe[16j and Asai[17j. They used a multivariate 
Student's t-distribution in order to cover the tails of the posterior density and determined the 
parameters to specify the distribution by using the ML technique. Here we also use a multivariate 
Student's t-distribution but determine the parameters through MCMC simulations. 
The (p-dimensional) multivariate Student's t-distribution is given by 
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where and M are column vectors, 
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and Mi = E[9i). S is the covariance matrix defined as 



= E[{e - M){e - Mf]. (17) 



V \s a, parameter to tune the shape of Student's t-distribution. When v ^ cc the Student's 
t-distribution goes to a Gaussian distribution. At small v Student's t-distribution has a fat-tail. 
Since eq. p5]) is independent of the previous value of 9, eq. p^ is used in the MH algorithm. 

To use eq. (|15p as a proposal density we have to know the values of M and S. We determine 
these unknown parameters M and S as follows. First we make a short run by an MCMC method 
and sample some data. Then we estimate M and S from those data. Second substituting the 
estimated M and S to eq. (jl5p we perform an MH simulation with the proposal density. After 
accumulating more data, we recalculate M and S, and update M and S of eq. p3]) . By doing 
this, we adaptively change the shape of eq. (fT5|) to fit the posterior density more accurately. We 
call eq. p3]) constructed in this way "adaptive proposal density". 

The random number generation for the multivariate Student's t-distribution can be done 
easily as follows. First we decompose the symmetric covariance matrix E by the Cholesky 
decomposition as S = LL*. Then substituting this result to eq. (fT5|) we obtain 



1 + 



(18) 



where X = L ^{9 — M). The random numbers X are given hy X = Y^J—, where Y follows 

N{0,I) and w is taken from the chi-square distribution u degrees of freedom Xu- Finally we 
obtain the random number 9 hy 9 = LX + M. 



7. Empirical Studies 

In this section we make empirical studies based on daily data of Nikkei 225, DAX and Hang 
Seng indexes. The sampling period is 4(2) (3) January 1995 to 30 December 2005 for the Nikkei 
225 (DAX) (Hang Seng) index. The index prices pi are transformed to returns as 

ri = WOln{p,/p^^l-s), (19) 

where s is the average value of ln(pj/pj_i). Fig[T] shows the time series of returns of the Nikkei 
225 index calculated by eq. (fT9|) as an example. 

The adaptive construction scheme is implemented as follows. First we make a short run by 
an MCMC method. Any MCMC method can be used. Here we use the Metropolis algorithm. 
We discard the first 5000 data as burn- in process (thermalization). Then we accumulate 1000 
data to estimate M and S. The estimated M and S are substituted to g{9) of eq. (fT5]) . The 
shape parameter u is set to 10. We re-start a run by the MH algorithm with the proposal density 
g{9). Every 1000 update we re-calculate M and T, using all accumulated data and update g{9) 
for the next run. We accumulate 100000 data for analysis. The results analyed are summarized 
in Table 1. 

We examine the efficiency of the algorithm by measuring correlations between sampled data. 
Figl2] shows the autocorrelation function (ACF) of the sampled a with the Nikkei 225 index 
data. The ACF is defined as 

. . i Ef-iixU) - {x)){x{j + t)- (x)) , ^ 

ACF(t) = ^^^"^^ -4r^ — —, 20 
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Figure 1. Time series of returns of the Nikkei 225 index. 
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Figure 2. Autocorrelation functions of a sampled with the Nikkei 225 index data. 



where (x) and are the average value and the variance of certain successive data x respectively. 
The ACF quickly decreases with t, which indicates that the correlation between the sampled 
data is very small. We also find the similar behavior for the other parameters. 

In order to analyze the autocorrelation quantitatively we measure the integrated 
autocorrelation time Tint given by 



(21) 



1=1 



"2Tint" is also called inefficiency factor. If no correlation exists in the data 2Tint takes one. The 
results of 2Tint are summarized in Table 1. We find that all the 2ri„t are very small, less than 2, 
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Figure 3. The diagonal element Vu as a function of Monte Carlo time. 



which indicates that the correlation between the data sampled by this method is small. 



Table 1. Results of QGARCH parameters. SD and SE stand for standard deviation and 
statistical error respectively. Theoretically there is a relation of SE ~ ^J^^^SD. Here statistical 
errors are estimated by the jackknife method. 





a 






7 


Nikkei 225 


0.07872 


0.89390 


0.06219 


-0.12403 


SD 


0.011 


0.013 


0.013 


0.021 


SE 


0.00003 


0.00004 


0.00005 


0.00007 




2.0 ±0.1 


2.0 ±0.1 


2.0 ±0.2 


1.8 ±0.1 


DAX 


0.09198 


0.89564 


0.03004 


-0.08483 


SD 


0.011 


0.011 


0.0064 


0.015 


SE 


0.00004 


0.00005 


0.00004 


0.00006 




1.77 ±0.06 


1.78 ±0.05 


1.80 ±0.07 


1.61 ±0.06 


Hang Seng 


0.07638 


0.91168 


0.03202 


-0.08678 


SD 


0.009 


0.0098 


0.007 


0.007 


SE 


0.00005 


0.00005 


0.00003 


0.00007 


2Tint 


1.80 ±0.06 


1.78 ±0.06 


1.79 ±0.06 


1.75 ±0.06 



Figl3] and m show the convergence property of the covariance matrix. Here F is a 4 x 4 matrix 
defined hy V = E[{9 — M){9 — M)*]. The elements Vu and V12 quickly converge to certain 
values. We also see the similar behavior for the other elements. 

FiglUshows the acceptance at the MH algorithm with the adaptive proposal density of eq. (ll5p . 
Each acceptance is calculated every 1000 updates and the calculation of the acceptance is based 
on the latest 1000 data. At the first stage of the simulation the acceptance is low. This is 
because M and S have not yet been calculated accurately at this stage as seen in figs. [SJIH The 
acceptance increases quickly as the simulation is proceeded and reaches plateaus of about 80%. 

The QGARCH model is capable of capturing leverage effects. In order to see the impact of 
leverage effects Pagan and Schwert[25], and Engle and Ng[26] proposed the use of the so-called 
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Figure 4. The off-diagonal element V12 as a function of Monte Carlo time. 
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Figure 5. Acceptance at the MH algorithm with the adaptive proposal density. 



news impact curve. The news impact curve is the functional relationship between conditional 
variance at time t and the shock at time t — 1, yt-i- 

The news impact curve of the QGARCH model is given by 



= uj + 7yt_i + ayt_i + /3cr^, 



(22) 



where is the unconditional variance given by 



15 



T 




- Nikkei 225 

- DAX 




-5 





yt-i 



5 



10 



Figure 6. News impact curves of the Nikkei 225, DAX and Hang Seng indexes. 

Figl6] shows the impact curves of Nikkei 225, DAX and Hang Seng indexes. Since the values 
of 7 for three indexes are all negative the three indexes exhibit the leverage effect. As seen in 
figiHthe news impact curve is higher for negative yt-i than for positive one. Since all the three 
news impact curves are very similar each other the three indexes turn out to have the similar 
responce to positive or negative news. 

8. Conclusions 

We have performed the Bayesian inference of the QGARCH model by the MCMC algorithm. 
The MCMC algorithm was implemented by the MH method with the adaptive proposal density. 
The adaptive proposal density is assumed to be the Student's t-distribution and the distribution 
parameters are determined by the data sampled by the MCMC simulation. The distribution 
parameters are updated during the MCMC simulation adaptively to match the posterior density 
of the model parameters. 

We have applied our method for Nikkei 225, DAX and Hang Seng indexes. We find that the 
autocorrelation times between the sampled data are very small, typically 2rj„f is less than 2. 
Thus our method is very efficient for generating uncorrelated Monte Carlo data. 

The QGARCH model is designed to capture the leverage effect. We find that all the three 
indexes show the leverage effect, i.e. the volatility is higher after negative observations than 
after positive ones. 

The acceptance of the MH algorithm quickly reaches a plateau of about 80% already at 
the beginning of the simulation. This means that the distribution parameters are evaluated 
accurately already at the beginning of the simulation, as shown in figsISHU This observation 
suggests that in practice one may stop the update of the parameters at some stage of the 
simulation. 
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